Brownian motion of free particles on curved surfaces 
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Brownian motion of free particles on curved surfaces is studied by means of the Langevin equation 
written in Riemann normal coordinates. In the diffusive regime we find the same physical behavior 
as the one described by the diffusion equation on curved manifolds [J. Stat. Mech. (2010) P08006]. 
Therefore, we use the latter in order to analytically study the whole diffusive dynamics in compact 
geometries, namely, the circle and the sphere. Our findings are corroborated by means of Brownian 
dynamics computer simulations based on a heuristic adaptation of the Ermak-McCammon algorithm 
to the Langevin equation along curves, as well as on the standard algorithm, but for particles per- 
pendicularly attached to the surface through highly stiff springs. The short-time diffusive dynamics 
is found to occur on the tangential plane. Besides, at long times and compact geometries, the mean- 
square displacement moves towards a saturation value given only by the geometrical properties of 
the surface. 

PACS numbers: 05.40.-a, 83.10.Mj, 82.70.-y 



I. INTRODUCTION 

Diffusive processes occur in a wide diversity of physi- 
cal areas; ranging from soft condensed matter to particle 
physics and astrophysics [TH3]. In the latter case, such 
processes are relevant within the framework of special 
and general relativity (see, e.g., Ref. [3j and references 
therein). Furthermore, in the last few decades, an in- 
tense activity has emerged in biophysics to study protein 
diffusion on curved surfaces [5]. 

The transport phenomena of proteins occurring inside 
cell membranes are interesting and complex since they 
determine the flux of nutrients between the cell and its 
exterior affecting, in consequence, the cell functionality 
[6j. From the theoretical point of view, it is difficult 
to describe the lateral diffusion of integral proteins or 
lipids mainly because the interactions with the remain- 
ing components of the membrane and the protein finite- 
size effects [3 |8] . In addition, there are membrane cur- 
vature contributions [5] and thermal fluctuations that 
produce shape undulations in the membrane and, thus, 
have a crucial contribution to the transport phenom- 
ena. Furthermore, the coupling between lateral motion 
of proteins and thermal shape fluctuations has been re- 
cently considered by Reister and Seifert jlOj , whereas pro- 
tein effects on the shape of the membrane surface were 
studied in references [HI [12]. Protein diffusion is also 
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affected when changes in the membrane thickness are 
present [T3J [T3]. The study of protein lateral diffusion 
has been made possible by assuming the membrane as a 
two-dimensional regular and curved sheet and projecting 
the Smoluchowski equation on such a surface [HI H5H2T] . 
In this approximation, both thermal shape fluctuations 
and finite-size effects have not been taken into account 
explicitly. 

Recently, one of us explored the curvature effects on 
the Brownian motion, at short times, of free particles 
when their movement takes place in a d-dimensional Rie- 
mannian manifold [22] . In that work, a general method 
is provided to derive all moments, (s™), of the probability 
density, P(x,x',t), of finding a particle at a point x and 
time t when it started its movement at x' . In particular, 
the method is shown in which the mean-square displace- 
ment (MSD) is deviated from the planar expression by 
curvature terms which are O (d)-invariant, as well as in- 
variant under general coordinate transformations. Here, 
s is the geodesic distance in the d-dimensional Rieman- 
nian manifold and (• • • ) stands for an ensemble average. 
It is also found that the second moment, up to second 
order in t, is given by |22j 

(s 2 ) = 2dD Q t - ^R g (D Q t) 2 + ■■■ , (1) 

where R g is the Ricci scalar curvature of the general Rie- 
mannian manifold and Dq is the free-particle diffusion 
coefficient. It is noteworthy to mention that the same 
curvature effects have also been explored to discuss the 
role of diffusion within the general relativity context [23J . 

Although the diffusion equation is suited to study the 
Brownian motion of free particles on curved surfaces 
[22] . a more complete description is provided by the 
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Langevin equation. The latter is based on the New- 
ton's equation of motion but including a rapidly fluc- 
tuating force, Gaussian distributed, representing the in- 
teraction among the particle and the solvent. It is well- 
known that in Euclidean open spaces the MSD calcu- 
lated from the Langevin equation reproduces the stan- 
dard Einstein kinematical relation. In this kind of spaces, 
both Langevin and diffusion equations describe the same 
dynamical behavior at the diffusive time regime, i.e., 
t ^> tb = M /£, where £ is the friction coefficient of the 
solvent, M is the particle mass and t b the momentum 
relaxation time 24 . In fact, it is found that using meth- 
ods of stochastic processes on manifolds, the diffusion 
equation can be derived from a Langevin-like equation 
[25j . Furthermore, Langevin- type equations have also 
been used to study Brownian motion in crystals with 
defects and including torsion effects [26], the lateral mo- 
tion of proteins [27] and diffusion in fluctuating ruffled 
membranes [2"8] . 

In this work, we write the Langevin equation for mani- 
folds without including torsion. Our starting point is the 
Newton's equation for free particles in a d-dimensional 
hypersurface. Free means here that particles do not in- 
teract between each other and non external force is acting 
on them, they are nevertheless restricted to move on the 
hypersurface. We find a set of equations that correspond 
to the global version of the equations found in reference 
[2"o] . The resulting equations allow us to study the cur- 
vature effects using the tools of the Riemann normal co- 
ordinates (29) . We reproduce the same leading curvature 
effects in the diffusion regime as in reference [22] , but the 
Langevin equation also allows us to study shorter time 
regimes. Taking all this together, we investigate the par- 
ticle dynamics for T so i vent -C t -C tb, for tb -C t < tq, 
as in reference [22], and for t 3> tq. Here, T so i vent is a 
characteristic time for the position and momentum re- 
laxation of the solvent molecules, at which the Langevin 
description is not longer valid, and tq is the time scale 
thereafter the effects of the geometrical properties of the 
manifold become evident. We explicitly analyze the dy- 
namics of particles confined along a circle, as well as on 
a sphere. 

We test our predictions by means of Brownian dynam- 
ics computer simulations based on a heuristic adapta- 
tion of the Ermak-McCammon algorithm to the Langevin 
equation along curves, as well as on the standard algo- 
rithm, but for particles perpendicularly attached to the 
surface through a spring-like force. In the first case, 
which is here only applied to the circle, the particles are 
allowed to move in any direction with equal probabil- 
ity, but the geodesic distances they travel are Gaussian 
randomly distributed. In the second case, the particles 
are allowed to move everywhere in the d + 1-dimensional 
Euclidean open space containing the hypersurface, but 
firmly attached to the latter by some kind of virtual 
springs acting perpendicularly to any point of the sur- 
face, at any time. In the limit case of very stiff springs, 
we get the same results from both routes, as we will see 



further below. 

The manuscript is organized as follows. In section II 
we present the Langevin equation for curved manifolds, 
written in both local and global coordinates. In addition, 
we study the curvature effects on the MSD at the follow- 
ing regimes: T so i vent < t < t b and t b < t < t g . In 
section III we study the particle dynamics on the geomet- 
rical regime (t> tg) by means of the diffusion equation 
on curved manifolds. In section IV we explicitly com- 
pare our predictions for particles restricted to move along 
a circle and on a sphere with Brownian dynamics com- 
puter simulations. Finally, in section V we summarize 
some concluding remarks and perspectives of our work. 



II. LANGEVIN EQUATION ON CURVED 
MANIFOLDS 



A. Global coordinates description 

We now specify the basis of our Langevin dynamics 
formalism. It is defined over an Euclidean hypersurface 
S C R d+1 , which is also defined as the points X e R d+1 
such that $ (X) = 0. The Langevin equation needs to 
include an holonomic constraint in order to bound a point 
particle to S. 

Let us denote the momentum P e 7x(5), where Tx(S') 
is the tangent space at the point X, and the position 
X £ R d+1 of the particle. From a classical point of view, 
the addition of the term A$(X) to the free-particle La- 
grangian allows us to impose an holonomic constraint on 
S. Indeed, the equation of motion is P = A V<I> (X) and 
the required constraint is $(X) = 0. We should remark 
that A = relaxes the constraint. For the Langevin 
equation defined on 5, we simply include the previous 
constraint, a friction term and a stochastic force 

P = -C P/M + A V$(X) +f (t) (2) 
X = P/M, (3) 
*(X) = 0. (4) 

The second term of equation ^ represents the force 
caused by the holonomic constraint. The stochastic 
force f (t) is chosen such that it satisfies the fluctuation- 
dissipation relations 



(/*(*)) = o, 

(fi(t)fj(t')) = ns i:j s(t-t'), 



(5) 



where (• • • ) stands for the average in the ensemble of 
field forces Gaussian distributed over R d+1 space (see Ap- 
pendix A). R d+1 is a copy of T x (5) xl, i.e., the ensemble 
is given by all possible configurations of forces belonging 
to K d+1 . In other words, for each point of the surface, 
we have a Gaussian distribution of forces. 

The Lagrange multiplier A can be obtained using the 
constraint (4]) as follows. A time derivative on this con- 
straint implies that 



V*(X)-P = 0, 



(6) 
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where V represents derivation in the space M <i+1 . Since 
the momentum P £ Tx(S'), then V$ is normal to the 
tangent space. Thus, the normal vector to the surface, 
i.e., normal to T X (S), is given by n = V$ (X) / |V$ (X)|. 
Second derivative on equation ([6| gives 



n P 



M 3 ' 



(7) 



where p a = d(g ab p b )/dt. The partial derivative d b e a 
can be computed using the Weingarten-Gauss equations 
V a e5 = — K ab n, where V a is the covariant derivative, 
compatible with the metric g ab , and K ab the component 
of the second fundamental form [30]. The Weingarten- 
Gauss equations can be written in terms of the Christoffcl 
symbols as d a e b = T c ba e c — K ba r\. By using this equation 
in the momentum time derivative we obtain 



with Gij = 9j9j$/|V$|. Now, we get A by equat- 
ing Q and the normal projection of (J2j). Then, A = 
-P i G ij P j /M |V$|-n-f/ |V$|. Therefore, the Langevin 
equation involves a non-linear term proportional to a sec- 
ond power in momenta, 



F 



c " ^K abP a p b n. (10) 



— Garvin = -— P 
M M 



The local coordinates version of the Langevin equation 
can be straightforwardly obtained by substituting equa- 
tion (10 1 in equation (Tsl) . Hence, the tangent projection 



P (f (t)) , (8) takes the form, 



where we have introduced a projector P that maps a vec- 
tor V e T X (S) x R = R d+1 into the tangent space. This 
projector is given by P y = <5 y — n % n> and it satisfies the 
usual relationship: P 2 = P. The quadratic term in mo- 
menta is not a surprise since the left-hand side of equation 
^ corresponds to the ordinary kinetic term for a par- 
ticle over a hypersurface. In other words, the Langevin 
equation ^ reduces to the geodesic equation when both 
the friction and the stochastic force vanish. This will be 
clear further below when we write equation ^ in local 
coordinates. 

We point out that the G matrix encodes the surface 
geometry. For instance, the constraint <&(X) = a • X + b 
defines a plane in Euclidean space, where a is a constant 
vector and b a real number. In this particular case, the 
G matrix is zero and the normal vector of the surface 
is constant, n = a/|a|, as it is required for a planar 
geometry. In the case of a sphere of radius R, we have 
$(X) = X 2 /i? 2 — 1 and the normal vector satisfies n = 
X/_R; the matrix G is given by Gij — Sij/R. 

It is important to mention that the Langevin equation 
written in global coordinates ^ may be helpful for nu- 
merical simulations. Additionally, the global description 
is also the natural starting point to introduce ambient 
interactions, where the extrinsic geometry may play a 
role. 



B. From a global to a local coordinates description 

We now provide a local coordinates description of the 
Langevin equation (jsj). In local coordinates a hypersur- 
face is parametrized by the mapping X : U C M 2 — >• S, 
where a particular point in S is given by p = X (£ a ), be- 
ing £ a the local coordinates (a = 1, • • • , d). In these coor- 
dinates, we have X = e a £ a , P = e a p a , and P (f) = e a f a , 
where p a = g ab Pb is the local momentum and e a = 9 a X 
the tangent vectors (note that d a — d/d£ a ). Thus, the 
first derivative of the momentum is given by 



1 

M 



d b e a p a p" + e a p a , 



(9) 



c 



P = ~M P ' M 
t - 9 ab Pb/M 



r c baP b P a + f 



while the normal projection is given by 



K v. — p 1 G p j 



(11) 



(12) 



Equations (11) and (12) are the local version of the 



Langevin equation ([8|) and they can be easily extended 
to the case of manifolds with torsion (see [26]). As we 
mentioned above, the quadratic contribution in momenta 
is just the geodesic contribution. The normal projection 
provides a geometrical identity that allows us to derive 
the extrinsic curvature in terms of the G matrix. This 
identity is not casual; it is actually the same found at 
the level set formulation of the hypersurfaces geometry 
[30] . The extrinsic properties of the hypersurface are not 
playing a crucial role, as it is expected from the very 
beginning, since the motion of the particle intrinsically 
occur on the hypersurface. 

Regarding the fluctuation-dissipation relations, the 
stochastic forces satisfy the following properties: 



(/-(*)) 

WW)) 



o, 

ns ab 5{t-t'), 



(13) 



where 5 ab is the two-dimensional Kronecker's delta. 
These relations are equivalent to their global version (see 
Appendix A). 



C. Dynamics beyond a local neighborhood 



Based on equation (11), it is clear that the particle 



dynamics does not depend on the extrinsic properties of 
the geometry. This means that the dynamics on a hyper- 
surface can be studied in a Ricmannian geometry, what 
we do in the following. Additionally, we mainly focus 
on the diffusion in the weak curvature regime. Let us 
recall that if V C S is a local neighborhood of S, the 
map X : U C M. d — > V is a local diffeomorphism |30j . 
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then V = X({7) = R d . This implies that in a local neigh- 
borhood, we should have the same particle dynamics as 
found in planar spaces (see reference [53] for the K 3 case) . 
Thus, it makes sense to study curvature effects around 
the Euclidean solution. 

In what follows, we first review the particle dynam- 
ics on the Euclidean geometry S — M. d , i.e., when the 
curvature is zero, and, second, we expand the Euclidean 
solution in order to study the leading curvature effects 
on the particle dynamics over the surface. 



1. Euclidean geometry S — M. d 

In the Euclidean geometry, both the global and local 
descriptions are the same; the Euclidean metric is simply 
9ab = S a b and the Chrystoffel symbols are zero. In this 
case, the Langevin dynamics formalism reduces to the 
well-known standard equations (53] 



P = 
t = 



M 



M 



P , 



(14) 



and their solution can be written as 

p c (t) = p c Q e-^ t + [ dt'f c (t')e 
Jo 

JL rt 

M 



v c {t) 



V6 



dt'p c {t'). 



(15) 



Proceeding along the same lines, we also get the MSD: 



<* 2 (t)> 



dSlM 



M 2 V 



2-4-* 



1) 



-2(1 -e"* 4 ) 



(19) 



In the diffusive regime, f > r^, the average kinetic 
energy reaches its equilibrium value. This allows the 
evaluation of f2 from the equipartition theorem. Thus, 
(p c {t)p c (t)) = dk B T/2 and Q = 2(k B T, where k B is 
the Boltzmann constant and T the absolute tempera- 
ture. We also observe that in this time regime the MSD 
reproduces the standard kinematical Einstein relation 
(s 2 (i)) — 2dD t, where D = k B T/( is the free-particle 
diffusion coefficient [21] ■ We should point out that the 
value of f2 is independent of whether the space is curved 
or not, since it only depends on quantities intrinsic to the 
fluid, as solvent friction and particle dimension. 

Higher order temporal correlation functions are 
also useful. In particular, we will see below 
that the four-points function G abcd (ti,t2,t 3 ,t 4 ) = 
(p a {ti)p b (t2)p c {t3)p d (t4)} is necessary in order to obtain 
the leading curvature corrections. This correlation func- 
tion can be computed by using the Wick's theorem [25] : 

G abcd {t x MMM) = (p a (h)p\t 2 ))(p c {h)p d {t A )) 
+ (P^h^ih)) ( P b (t 3 )p d (t 4 )) 

+ (p a (h) P d (t 4 )) ( P b (t 2 )p c (t 3 )po) 

2. Leading weak curvature effects 



Averaging equations (15 1 over the ensemble of stochastic 
forces, one easily obtains 



(y c (t)) = y c + ^o(i-e-**) 



(16) 



We observe that the mean momentum decreases expo- 
nentially with time (with the decaying time scale t b — 
M/Q and the particle position is shifted by p c (0)/<^ at 
long-times. 

We now consider for simplicity that = and = 0. 
Other physical quantities of interest are both the mean 
quadratic momentum, i.e., (p c (t)p c (t)) , and the mean 
square displacement, s 2 = y c y c - In order to calculate 
them, it is useful to find the temporal correlation func- 
tion between two momenta, p a {t) and p b (t'), at times t 
and t' , given by (see Appendix B for further details) 



{P a {t)p\t')) 



VIS 



ab 



.-£lt-t'l 



-(*+*') 



(17) 



Using previous equation, it is straightforward to obtain 
the mean quadratic momentum: 



(p C (t)Pc(t)) 



dClM 



(l-e" 2 *') 



(18) 



We turn now to the derivation of the leading weak cur- 
vature effects on the particles dynamics. As we already 
discussed, the Langevin equation is quadratic in the mo- 
mentum and that contribution is coupled to the particles 
positions through the Chrystoffel symbols. The result- 
ing equations are difficult to solve analytically, among 
other reasons because the left-hand side of equation ( 11 ) 
involves a temporal derivative of the metric. Using 
p a = g ab pb (p a is independent of the metric), the lo- 
cal Langevin equation allows us to obtain the following 
expressions: 



Pd 



y" 



1 
M 



Pd 



1 

J] 



fled 



t^c bf ah , p 

^3cdl ba 9 9 PfPh + Id, 



M 



ab ^ 

9 Pb- 



(21) 



In order to explore curvature effects, we expand equa- 



tion (21| around the planar solution (15 1. To reach this 
goal, we use the Riemann normal coordinates |29| . In 
normal coordinates, we have 



9 



a b 



S ab 
1 



l -R a cd b y c y d + 0(y 3 ), 



hi I a 



R c adb )y d + 0(y 3 ), (22) 
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where R a bcd are the components of the Riemann curva- 
ture tensor. Using (22) in (21) one obtains 



Pd 



y" 



M 



AI 



1 , , 

Pd - T^RdbfaVh < P a + ■ 
/rafc ^ r>a bed, \ 

- ~R cd y v h — )Pb 



fd 



(23) 



We should notice that the Langevin equation in Eu- 
clidean geometries ( 14 ) is recovered when the curvature 



vanishes. In order to find a solution around the Euclidean 



case ( 15 ) we expand the momentum and position in the 
following way: Pd — qd + 8qd and y a = z a + Sz a , where 
qd and z a are the solutions for zero curvature, given by 
(15). Here, we have assumed that Sq = and 5z = 
when R a hr .j = 0. If we consider only linear terms in cur- 



(24) 



bed 

vature, we obtain the equation for Sqd'. 

C 1 f b 

S 1d = ~M 6 1d - 7^jRdbf a z J q q a - 

The second term of the right-hand side does not depend 
on 5q; it depends only on time. The integration of equa- 
tion ( 24 ) is similar to the one in the planar case. The ini- 
tial condition for 5qd{t) is Sqa(0) = 0, since pd(t) satisfies 
Pd(fi) = 9d(0). Therefore, the momentum, up to linear 
terms, in an arbitrary Riemannian geometry is given by 



Pd(t) = qd(t) 



3M 2 



Rdbc 



dt't 



t" 



< / dt"q c {t")q b {t')q a {t% 
>0 



(25) 



and the position, up to linear terms as well, takes the 
form 



y a (t) = z a {t) 



i 

3M 3 ~ 



t" 



R a bcd / dt' / dt"e-i(*'-*") 



t" 



x / dt"'q b {t")q c {t"')q d {t"). 



(26) 



Wick's theorem allows us to determine the temporal cor- 
relation functions of q a (t). Therefore, we have found that 
the odd correlations vanish, as in the case of the mean 
values of the momentum and position: (p a (t)) = and 
(lla CO) = 0- This means there is not preferential points 
on the surface and the mean values are independent of 
the geometry. This result may change however for non- 
zero initial conditions. 

Up to linear terms in the curvature we obtain the fol- 
lowing expectation value for s 2 (t): 



The quantity J a ijk(t) captures the dynamical contribu- 
tion that appears in the weak curvature regime. In addi- 
tion, the four-points correlation function G ai jk is defined 
according to equation ( 20 ) ; it is built by the products of 



two-point correlation functions and each of them carries 
a Kroncckcr's delta. Hence, using the symmetries of the 
Riemann tensor, the MSD reduces to 



( s *(t)) = (z a (t)z a (t))-^J(t) + 



(29) 



where R g is the scalar curvature and 



J(t) = dhj dt 2 J dt 3 J dfce"* <*»"*») 

x {G(h,U)G(t 3 ,t 3 ) - G{hM)G{UM)) ■ (30) 



*3 



Equation (30) can be straightforwardly integrated (see 
Appendix B). Equation p9| represents the MSD 
(geodesic mean square displacement) in the weak cur- 
vature regime. As we can appreciate from equation ( |30[ ) , 
the time scale tb = M/£ defines two time regimes: The 
one with i < rg (but very much larger than T so i ven t) 
or the ballistic regime, and the one with f > or the 
diffusive regime. In the first case, the MSD is given by 



<* 2 (t)> 



2d(k B T 
M 2 



(31) 



The cubic term is the ordinary contribution to the ballis- 
tic regime when the initial condition is p$ = (it becomes 
quadratic in t for non-zero initial conditions [24 ). The 
next curvature contribution is of order i 6 ; typically neg- 
ligible unless there is a region of very high curvature. 

In the diffusive regime, t > rg, the function J (t) re- 
duces to J (t) pa (Dot) 2 . Therefore, the MSD becomes 



(s\t)) = 2dD Q t- 2F ^{D tf 



(32) 



This result is the same found in reference |22j by means 
of the diffusion equation on curved manifolds. The MSD 
shows a deviation from the planar result due to curva- 
ture effects. Furthermore, equation (32) also shows the 



(s 2 (t)) = (z a (t)z a (t)) 



where (z a (t)z a (t)} is the same as in equation (19), 



Jaijk(t) = dti dt 2 



dU 



raise of two different diffusive regimes: The one with 
tb t < tg, and the so-called long-time regime, also 
called geometric regime, t^>r G . Here, r G = 3d/ \ R g \ D 
stands for the time thereafter the curvature effects be- 
come dominant. Thus, up to linear order in the cur- 
vature, the Langevin equation and the diffusion equa- 
tion are equivalent in the short-time diffusive regime. 
Nonetheless, additional curvature effects may be encoded 
in the Langevin equations for i > Tg, which we do not 
discern in this paper. 

It is important to mention that in the planar case, 
i.e., |i2 g | — > 0, the particle cannot feel any effect asso- 
ciated with the geometry (tq is never reached, then it 
dt/ueT w(' 2- * 3 ) grows towards infinity). Additionally, we should empha- 

size that in the particular case of d = 1 the MSD may 
(28) exhibit deviations that cannot be associated to R g , since 



,(27) 



() 



the Gaussian curvature of lines is zero. In fact, as we will 
see below, those effects are associated with the finite-size 
of the phase space. 

From now on, we will consider that the complete diffu- 
sive regime is well described by both the Langevin equa- 
tion and the diffusion equation. In the following section 
we discuss some properties of the diffusive motion of the 
particles along a circle, S 1 , and on a sphere, S 2 . 



with "f* being the complex conjugate of ^. We note that 
degeneracy of eigenvalues is explicitly considered in the 
sum. 

Now, let us consider an arbitrary observable O. Its 
dynamical behavior can be obtained using the formal ex- 
pression for P (x, x',t). The expectation value (O (x)) 
has a generic form; its structure around the geometric 
regime is determined by the smallest eigenvalues. Then, 
it can be written as follows: 



III. DIFFUSION IN S 1 AND S 2 

We now choose the diffusion equation in order to study 
the geometric regime (t tq) in the manifolds S 1 and 
S 2 (for a discussion on the diffusion on arbitrary hyper- 
spheres see, for example, reference |3I|). The diffusion 
equation on curved manifolds can be written as [22], 



(O (x)) ~B|)- a\e 



-DotXi 



(37) 



dP (x,x',t) 



D A g P{x,x',t), 



P(x,x',0) = — (5 (d) (x - x') , (33) 

where P(x 1 x',t)dv is the probability of finding the dif- 
fusing particles in the volume element dv — y/gd d x, given 
that they began to move at x' . The probability density 
distribution P (x,x',t) is normalized with respect to the 
volume v of the manifold and Dq is the free-particle dif- 
fusion coefficient. The operator A g , called the Laplace- 
Beltrami operator, is defined by 



f 



A g f = -^d a (Vgg ab d b f) 



(34) 



with g — det (g a b) and / is a scalar function. The ge- 
ometry is coupled to the Brownian motion through the 
metric g. The diffusion equation (33 1 is the same as the 



heat kernel equation and it has a lot of applications in 
the context of field theories on curved spaces |32j . 

The expectation value of a scalar function O defined 
on the manifold is given in the standard fashion, i.e., 



(0(x)) = / dv 0(x)P(x,x',t), 



(35) 



and (O (x)) depends on the initial point x' . The char- 
acteristics of observables in manifolds are related with 
the particular structure of P (x, x',t). Furthermore, the 
probability density distribution P(x,x',t) can be deter- 
mined by solving the eigenvalue problem — Agty = £ , <J , 
where E is the eigenvalue corresponding to the eigen- 
function ^ . In addition, it is known that for compact 
manifolds, the spectra of A g is discrete and it can be 
written in a growing sequence {Ao = 0, Ai, A2, . . . }, where 
A/ + i > A/ [33]. We also have a sequence of orthogo- 
nal eigenfunction ^1, #2, • • • in L 2 (S) (square-integrated 
functions of S). In this sense, the probability density dis- 
tribution can be formally written as |34| 

P (x, x', t) = J2 e- XlDot ^} {x') (x) , (36) 
1 



where a — A J dvO and aj = A J dv 

We can easily obtain some properties of any observable 



by looking at the particular form of equation ( 37 1 . For 



example, at long times it converges to the expectation 
value ao. This is a consequence of the compactness of 
the manifold. In physical terms, the positions of the par- 
ticles relax and their distribution does not longer evolve 
with time. The expectation value (O (x)) becomes a as 
a consequence of the finite size of the space. The quantity 
ao is the geometrical average of O; this is the reason we 
called this regime the geometric regime. Although coun- 
terintuitive, the values of the observables do not depend 
on the temperature for t 3> tq; it is only a function of 
the surface geometry. 



A. Brownian motion over S 

Brownian motion on the circle represents, after the 
motion on the straight line, the simplest example where 
there is a clear manifestation of the geometrical effects on 
the particle dynamics, but it is also the most fundamental 
one, since it is fully described by a single physical vari- 
able. It is also relevant for the theoretical and experimen- 
tal study of single-file diffusion in quasi-one-dimensional 
interacting systems (see, e.g., [35] and references therein). 

The circle is the mappping X : [0, 2ir] — > M 2 , where 
X = (i?cos tp, Rsmip), with R being the circle radius. 
The Laplace-Beltrami operator in this case takes the form 
A51 = -jp ^2 • The eigenf unctions of this operator form 
the complete orthonormal set {e' mf | m € Z} in L 2 (5' 1 ) 
and their corresponding eigenvalues are A m = —m 2 /R 2 . 

In order to study Brownian motion on S , we choose 
the following initial and boundary conditions: </?'(0) = 
and P (<p, 0, 0) = S (ip) /2itR. After some simplifications, 
the explicit solution of the diffusion equation is 



2-kR 



2E 



l2 i¥ cos (imp) ) . (38) 



In this case, the distribution is normalized with the 
perimeter of the circle, i.e., Jjds P(tp,t) = I, where 
/ = (— 7r, 7r) and ds = R dip. The distribution is also 
symmetric under the interchange ip — > —tp. 

The first moment, (s(t)), and the second moment or 
MSD, (s 2 (t)), of the distribution can be straightfor- 
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motion on the sphere in the framework of the diffusion 




tD /R 



FIG. 1: Mean square angular displacement as a function of 
time for free Brownian particles diffusing along a circle. The 
line corresponds to our theoretical result given by equation 
( 39 1 , and the symbols to the Brownian computer simulations 
results obtained by means of both the standard Ermak and 
McCammon algorithm (circles) and its heuristic adaptation 
to curves (triangles) . The error bars of the simulation data are 
smaller than the size of the symbols. There is no appreciable 
difference between the results. The straight lines stand for 
the short and long-time limits. 



wardly evaluated. The former is zero, since the distri- 
bution is an even function, whereas the MSD has the 
form 



(s 2 (t)) 
R 2 



OO 

m—l 



(-i) r 



(39) 



with s = Rip being the ar c-le ngth. On the one hand, the 
MSD given by equation Ety reduces to (s 2 (t)) = 2D t 
for short times (tb <C t <C Tq). On the other hand, at 
long times (t > t G = R 2 /D ) we have (s 2 (t)) = tt 2 R 2 /3. 
In the geometric regime the dependence is only on the 
size of the circle. The numerical evaluation of equation 



(391 is shown in figure I. 

As we mentioned previously, although the MSD in ( 39 ) 



deviates from the planar result, this difference is due to 
the finite size of the circle and not to curvature effects. 



We compare the predictions of equation ( 39 ) with com- 
puter simulations in figure 1. 



B. Brownian motion over S 

Brownian dynamics on the sphere has been studied 
by several authors using different approaches [TSl \T§\ [211 
|3BJI!IZ|- However, our choose to define displacement as 
the geodesic distance between two points on the surface 
provides complementary information. Nevertheless, wc 
do not discuss here the differences between alternative 
approximations. In contrast, we focus on the Brownian 



equation (33) 



The geometry of a sphere is encoded into the metric 
given by 



ds 2 



g ab dx a dx b = R 2 



sin 2 9dip 2 



(40) 



where R, 9 and tp are the radius, polar and azimuthal 
coordinates of the sphere, respectively. The Laplace- 
Beltrami operator on the sphere has eigenvalues and 
eigenvectors given by = £ (£ + I) and {Y( m (9, ip)} with 
I = 0, oo and m — —I, ■■■,£; Yi m (9, tp) being the stan- 
dard spherical harmonics. 

We choose x' to be on the north pole and take advan- 
tage of the rotational invariance. Besides, the boundary 
condition (33) is explicitly taken into account. The solu- 



tion of the diffusion equation is then 

°° of + 1 
P (^)=E^^ P/(coB0)exp 



£=0 



D t {I H 
R? 



1), 



(41) 



where Pi is the Legendre polynomial of order I. As in 
the previous case, we look for the information provided 
by (s(t)} and (s 2 (i)), but we have now s = R9. 

By means of the operator method defined in [22] it is 
possible to show that the short-time behavior of the MSD 
is given by equation ( 32 ) with the Gaussian curvature of 
the sphere, R g = 2//P7lt is interesting to note that the 
terms in the MSD that depend on the Gaussian curvature 
are always negative. This means that curvature effects 
only contribute to reduce the particle diffusion with time. 

In the geometric regime, t ^> tq = 3R 2 /D , we get 
from equation (37) 



R 



3 , D o* 
1 - - e ~ 2 ~^ 



R 2 



-4 



I - 



3tt 2 
4tt 2 - 16 ' 



(42) 



At the beginning the particles move around their ini- 
tial position, i.e., the north pole. After a long time, 
very much larger than tq, the expectation values (s) 
and (s 2 \ move towards the saturation values ttR/2 and 
(it 2 — 4)i? 2 /2, respectively. The particle have visited all 
the points on the surface and confinement dominates en- 
tirely the diffusive behavior; the saturation values only 
depend on the size of the sphere. The numerical evalua- 
tion of equation (42) is shown in figure 2. 

The expectation value of any observable O — 0(9,(p) 
on the sphere can be written as 



<O(0> ¥>)> = ? £(2*+l)<to(*)e- 



D e(i+i)t/R 2 



(43) 



e=o 



where go is the projection of O along the basis of Legen- 
dre polynomial. We explicitly show the functional form 
of go in Appendix C, for both g s and g s i. 




tD /R 



FIG. 2: Mean square polar angular displacement as a func- 
tion of time for free Brownian particles diffusing over a sphere. 
The line correspond to our theoretical result given by equa- 
tion (421, and the symbols to the Brownian computer simu- 



lations results obtained by means of the standard Ermak and 
McCammon algorithm (circles). The error bars of the simu- 
lation data are smaller than the size of the symbols. There 
is no appreciable difference between the results. The straight 
lines stand for the short and long-time limits. 



value zero and a covariance matrix given by the elements 
{SXniSXpj) — 2D°.£.At; these are the requirements 
needed to satisfy the fluctuation-dissipation theorem ([5| . 
The indices a and (3 run over the particles, and the in- 
dices i and j over the cartesian coordinates. In our case, 
we do not consider HI and, therefore, D a .a. = 5ij8 a p 



where Dq is again the free-particle diffusion coefficient. 
With this assumption, the second term in the right-hand 
side of equation ( 44 ) disappears and allows us to simplify 



drastically the calculation of the third and fourth terms 
of the same side. 

As we mentioned previously, the algorithm of Ermak 
and McCammon describes the temporal evolution of Eu- 
clidean variables. However, it can be still used to de- 
scribe the dynamics of particles on curved surfaces. This 
can be done by considering an external force that con- 
strains the movement of the particles on the surface. We 
demand that such a force does not contribute to the tan- 
gent displacements of the particles, i.e., this force has to 
act normal to any point of the manifold at any time to 
guarantee that it does not perform work on the system. 
Then, the simplest force that satisfies such requirements 
can be written as 



-fc(|X Q |-i?)n Q 



(45) 



IV. BROWNIAN DYNAMICS SIMULATIONS 
ON CURVED SURFACES 

A. Standard Ermak and McCammon algorithm 

In 1978, Ermak and McCammon introduced a method 
for simulating the Brownian dynamics of particles [55] , 
This method, which has been adapted in Euclidean co- 
ordinates, was derived from the Langevin equation and 
became consistent with the Fokker-Planck equation. Fur- 
thermore, such a method can be straightforwardly ap- 
plied when either hydrodynamic interactions are consid- 
ered explicitly or external forces act on the particles. 
This method has been successfully employed to study 
the structural and dynamic properties of a large variety 
of complex fluids, i.e., colloids, polymers, etc. [3H] 

The algorithm of Ermak and McCammon j38j is given 

by 



E U1J al3 



0=1 



N 
(8=1 



pF^At + SXa, (44) 



where N is the number of particles, (3 — (fc^T) -1 is 
the inverse of the thermal energy, ks being the Boltz- 
mann constant and T the absolute temperature. The 
hydrodynamic interactions (HI) are included through the 
diffusion tensor Dj^g, Fjg is the total force exerted on 
the /3-th particle and the index tells us that the vari- 
able must be calculated at the beginning in time at ev- 
ery step. The term 5~K a represents a random displace- 
ment with a Gaussian distribution function with mean 



where k is a spring-like constant, whose value is chosen 
in such a way that the particle displacements in the per- 
pendicular direction to the surface is basically negligible, 
R is the radius of either the circle or the sphere and n a 
is a unit normal vector. Hence, equation (45) is incorpo- 



rated in the standard algorithm for Brownian dynamics 
described in equation ( 44 ) to analyze the diffusion on the 
given surface. 



We should mention that the addition of force ( 45 1 into 



equation ( 44 ) has the same effect on the particle dynam- 



ics as the second term of the left-hand side in equation 
([8]), i.e., it only constrains the motion of the particles on 
the manifold. Thus, this kind of trick allows us to study 
the diffusion on curved surfaces through the use of the 
standard Ermak and McCammon algorithm. 



B. 



Heuristic adaptation of the Ermak and 
McCammon algorithm to curves 



Equation ( 44 ) takes the simple form X Q = X° + 5X Q 
for free particles, with (5X ai 5X aj ) = 25ijDoAt. In a 
d-dimensional Euclidian open space this process is equiv- 
alent to allow the particles to move in any direction 
with equal probability, as long as the distances they 
travel are Gaussian randomly distributed with variance 
(Ss a ds a ) = 2dD At. This is however the short-time 



behavior of the MSD in d-dimensional manifolds ( 32 1 



Hence, we hcuristically extend the Ermak and McCam- 
mon algorithm to curved manifolds by allowing the parti- 
cles to move in any direction with equal probability, but 
the geodesic distances they travel are Gaussian randomly 
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distributed, i.e., s = s° + 5s with (SsSs) — 2dDoAt, as 
long as t b < At < r G . 

In the particular case of a circle, this idea leads to the 
following algorithm: A uniform random number is gen- 
erated in the interval [0, 1]; the particle in turn is allowed 
to move in the clock-wise direction if the result falls in 
[0,0.5], otherwise the particle moves in the opposite di- 
rection; a Gaussian randomly distributed number with 
variance (Ssds) — 2DoAt is then generated in order to 
determine the arc-length the particle travels; these steps 
are repeated for every particle, many times, in order to 
construct the dynamics of the system in its natural se- 
quence. In our simulations, we let 1000 free particles to 
move in very short time steps, until they approximately 
cover a distance of 100 times the perimeter of the cir- 
cle. The large number of particles allows to improve the 
numerical precision of our results. 

We expect, on the one hand, the short-time behavior 
(As 2 (i)) = 2Dot, since this is included in the construc- 
tion of the algorithm. On the other hand, for very long 
times (t ^> tg) the particles has to distribute uniformly 
along the perimeter of the circle. Therefore, the geomet- 
ric behavior of the MSD must be given by the simple 
average of the geodesic square displacement 



(As 2 ft»r G )) 
R 2 



1 

2^ 



2n 



(^-M) 2 ^= y , (46) 



which agrees with equation (39). These and the interme- 
diate values of (As 2 (i)) are shown in figure 1. 

The extension of these ideas to the general case of 
curved surfaces will be presented elsewhere. 



V. CONCLUDING REMARKS AND 
PERSPECTIVES 

In this work the diffusion of free particles on curved 
surfaces is studied. After writing the Langevin equation 
and the fluctuation-dissipation theorem for curved sur- 
faces, we solved the former in the Riemann normal coor- 
dinates for weak curvatures, i.e., up to linear terms in the 
Riemann curvature tensor. From this solution, the dy- 
namics of the particles can be clearly separated in three 
regimes; the ballistic one, T so i vent <C t <C tb, and two 
diffusive regimes; short times, tq <C t < tq, and long 
times, or geometric regime, t tq. No effects of the 
geometry were found, neither in the ballistic nor in the 
short-time diffusive regimes. We therefore conclude that 
the local dynamics occurs in the plane tangent to the 
surface. Nevertheless, in the long-time diffusive regime 
only the geometric effects take place. The free particle 
diffusion coefficient Dq might be understood in terms of 
the short-time limit of the mean geodesic square displace- 
ment, (As 2 (tb <C t <C tg)) = 2dD t, in a similar way as 
in the case of interacting particles. The geometry then 
appears as an external force acting on the diffusing par- 
ticles, which can be recognized in the second term of the 
left side of equation d8|. 



We should remark that in the short-time diffusive 
regime the Langevin equation was found to have the same 
solution as the diffusion equation on curved surfaces [32] . 
We therefore used the latter in order to study the whole 
diffusive dynamics, although we do not discard the possi- 
bility to find differences between both equations beyond 
the weak curvature approximation. We particularly stud- 
ied the diffusion of free particles along a circle, S 1 , and 
on a sphere, S 2 . We do not expect curvature effects in 
S 1 since its Gaussian curvature R g is zero. However, the 
MSD displays a geometric diffusive regime due to con- 
finement effects, since the particles are unable to move 
beyond the region were the circle is placed. In S 2 the con- 
finement and curvature effects act together to define the 
geometrical regime. The difference between curvature 
and confinement effects is subtle and somehow counter- 
intuitive. This will be carefully reported somewhere else. 

In this work we also report some results from Brown- 
ian dynamics computer simulations. We obtained them 
by implementing the standard Ermak-McCammon algo- 
rithm, as well as its heuristic adaptation to curves. In the 
first case, we assumed that the particles are attached to 
the surface through springs acting on the particles every- 
where, always perpendicularly to the surface. The spring 
constant was adjusted to guarantee the particle dynam- 
ics very close to the surface. In the second case, which 
was only applied to the circle, we allowed the particles 
to move in every direction along the curve, every time 
displacing geodesic lengths given by random Gaussian 
number with variance (SsSs) = 2dD Q At, where d = 1 is 
the dimension of the circle. The quantitative compari- 
son of the theoretical results with the simulation data is 
shown in figures 1 and 2. 

A natural extension of this work could be the study of 
interacting particles. However, the interaction may pro- 
duce colored distributions for the stochastic forces in the 
Langevin equation [40 , so that Wick's theorem, which 
is of central importance in our calculations, were not 
longer valid. Nevertheless, it could be longer applied as 
an approximation, in the sense that the n-time correla- 
tion functions may be decomposed in terms of two-time 
correlation functions. In addition, both implementations 
of the Ermak and McCammon algorithm may be further 
used for interacting particles, as well as for other physical 
circumstances. 

For instance, in the case of the sphere, let us think in 
a cylindrical Brownian particle located at its centre. Let 
us assume that the particle is fixed to that point, but free 
to rotate around it. The particle only feels the solvent 
in which it is suspended, i.e., the motion of the particle 
can be described by the rotational Langevin equation for 
free particles |41j . If we imagine a virtual straight line, 
coincident with the longitudinal axis of the cylinder, but 
extending its length indefinitely in both directions, then 
the points where this line cuts the surface of the sphere 
may display a motion equal to the one we are generat- 
ing with our Brownian dynamics algorithms. Thus, the 
problem of the diffusion of free particles on a spherical 
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surface may be equivalent to the problem of the rota- 
tional diffusion of free particles. This equivalence could 
be extended to other surfaces, the circle for instance, as 
well as to the case of interacting particles. The latter 
including torques, which should be compatible with the 
interaction, given all particles are fixed to the centre of 
the sphere. No doubt that the study of this analogy is 
highly interesting. 



Appendix A: Fluctuation-dissipation theorem 



The stochastic force is Gaussian distributed for each 
point on the surface S. In global coordinates this dis- 
tribution is given by [25] 



d» = [I [dfi] exp {--L Jdt <)"/,/.. | 



(Al) 



where [dfi] is an appropiate functional measure. This is 
equivalent to a Gaussian vector field theory in one di- 
mension. The expectation values are defined by (O) = 
J dfi O I [dfi. In particular, the fluctuation-dissipation 
theorem |5]) can be verified using (Al). 



The force distribution (Al) also determines the 



fluctuaction-dissipation theorem in local coordinates 
(13). To show this, let us separate the force in tangent 
and normal components. Since, e a and n are given for 
each point, thus f = e a f a + n/ n is a biyective trans- 
formation between {/j}, with i = 1,2,3, and {f a ,fn}, 
with a — 1,2. Thus the measure ]X =1 [dfi] transforms 

to ri a =i \dfa] [df n ] J, where J = n- (ex x e 2 ) is the Jaco- 
bian. In addition, the argument of the Boltzmann weight 
can be splitted in these coordinates. Then, the measure 
dfi can be written as 



dfi = Y[ Wa] [dfn] J exp 



a=l 



29. 



dt{g ab f a f b + fl) 



(A2) 



Now, since the hypersurface is locally a plane we can 
always choose e a such that g a {, — S a b- Therefore, 



the local fluctuation-dissipation relations ( 13 ) can be 



straightforwardly obtained from (Al) 
detail allows us to establish that 



This technical 



30th global and local 
versions of the Langevin equation on curved surfaces are 
equivalent. 



Appendix B: Correlation functions 

1. Green function 

The correlation of two momenta for zero initial conditions 
can be computed from 



(v a (t) P b (t')) = 



dti / dt 2 e 
lo Jo 

x <r(*i)/ 6 (*2)>. 



(Bl) 



Next, we use the fluctuation-dissipation theorem (|13j). 
Thus the integration over variable t 2 leads to the follow- 
ing result 



f dt 2 e t2 l TB 5{t 2 -t 1 ) = 6{t l -ti)e 41 
Jo 



/tb 



(B2) 



where 8(x) is the Heaviside step-function. The remaining 
integral over t\ can be done for two cases t' > t and 
f < t. If t' > t then t' > t t for all t x e [0,i], therefore 
Bit' — t\) = 1. Now, if If < t then the integration for t\ 
can be splitted in two parts 

f dt x 6 (*' - tx) e 2tl / Tfl = f dtiQ (*' - t x ) e 2tl ' TB 
Jo Jo 

+ f dhOit 1 -h)e 2tl/TB . 



(B3) 

In the first integral t' > t±, since t\ £ [0, t']. Then for 
this integral it' — tx) = 1. For the second integral, we 
have t' < tx, since tx G [f , t]. Therefore 9 (f - tx) = 0. 
Now, joining these results and performing the elementary 
integrals we reproduce equation (17). 



2. Calculation of J function 

The determination of J{t) can be obtained from the 
calculation of 



J(t) 



1 



t rt rt 2 

dtx / dt 2 I dtp 



dUe~^~ H) 



x (G(tx, U)G(t 3 , h) - G(tx,t 3 )G(t 4l t 3 )) 



(B4) 



We should remark that the integral / (t, t') = 
f drG(r, t') appears in various places in the multiple 
integral ( |B4[ ). Thus, the function (B4) can be written as 
follows 



J(t) 



1 



pt2 

, , dtn 

M*J 2 J 



dUe~^ {t2 ~ t3) 



x ( G(t 3 , h) J dUl(t, U) ~ /(*, t 3 )I(t 3 , t 3 ) 

(B5) 



11 



The advantage to write J(t) in terms of I(t,r') is that 
t' < t. For the calculation of the function I(t, t') it is 
convenient to use the following equivalent expression for 
the Green function 



g(m') = r B n 



e -B0(t- t') sinh 



TB 



+ e t b0 (<' - t) sinh 



(B6) 



Performing its integral we obtain 



1 - 2e- 



(B7) 



No w, w e carry out the elementary integrations involved 
in ( B5 ) . We then get the following expression 



Jit) 



(r s D) 2 {e *a [, 



8 + e" 



t if t 

i— +e-s -37+12 — 



8 + 15-40e-E 



8e tb +36e t b 



+ 



+ 6^ (±-A 

TB V TB 



+ e 

t 



+ e 



TB 



1 + e t b 



(32 



(B8) 



Appendix C: Expectation values for Brownian 
motion over S 2 



The expectation values for O — O(0, <p) can be calculated 
from 



(o(e^)) = -J2m + ^)9o (£)> 

where 



-Dl(l+l)t/R 2 



(CI) 



£=0 



f 7T /■27T 

(£)= I / dflc^ sin 6» O (6, ip) P e (cos 9) . (C2) 



Equation (C2| depends explicitly on the chosen form of 



O. In general, equation (CI) cannot be written in a 



closed form and it must be studied numerically. In par- 
ticular, we discuss here the mean values of the functions 



O 



Re, and O 



In order to have a more 



manageable form for these expectation values we use the 
following identity 



p,( COS 0) = (-r/£ 



fc=0 



e J k \cos[(£-2k)6}. 



(C3) 

Now, in order to obtain (C2) we perform the integration 



for even and odd vaues of I. After performing the elemen- 
tary integrations, we obtain the following results. For 
s = R9, g s (£) is zero for even values of £, and for odd 
values of I it takes the form 



5s (2p+l) = 



_ l 

2 

P 



_ 1 

2 

P+l 



2p+l 

n J? (2(p-fc) + l) 2 



2p + 1 - k 
1 



(C4) 



where the last sum does not take the values k = p and 
k = p + l. For s 2 = R 2 8 2 , it is not difficult to show the 
identity g s 2 (2p + 1) = irg s (2p + 1) for odd values of I. 
However, for even values of £ we find 



2p 



g s 2 (2p) = 



fe=0 



2 P+ Lk^'- 



(C5) 



where H is a function defined as 



H(z) 



_ 12z 2 +4-tt 2 (z 2 -iy 



and 



= x (x - 1) (x - 2) ■ ■ • (x - n + 1) jn\ 
is the binomial coefficient [4"2] . 
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